%to read in edi files
% compute the MT parameters
%conour / interpolate them and display as a map

%latest date 8.5.3

[file,path] = uigetfile('*.edi', 'Directory for EDI Files');

K=dir(path);

nfiles = length(K);

for i = 3:nfiles,
       fp=fopen([path K(i).name],'rt');
       if fp==-1,
          return;
       end;
       [SPM,er]=edi_in(fp);
       
       bd4 = find(SPM.head.frq<0.1875); % band 4
       if length(bd4)>0,
          frq=SPM.head.frq(bd4);
          for j = 1:length(bd4),
             SPMat(j,:,:)=SPM.spectra(bd4(j)).data;
          end;
          Lat(i)=SPM.head.Lat;
          Long(i)=SPM.head.Log;
          C1 = pcoh(SPMat,'Ex');
          C2 = pcoh(SPMat,'Ey');
          CSite(i)=sum(C1 + C2)/(2*length(bd4));
       else,
          CSite(i)=0;
          fprintf('No data for site %s\n',SPM.head.site);
       end;
       
          SS(i,:)=sprintf('%5.3f',CSite(i));
                
          clear C1 SPMat;
  end;
  Long(1:2)=[];
  Lat(1:2)=[];
KK= find(Lat>9&Lat<13);
 plot(Long(KK),Lat(KK),'k.');
 text(Long(KK),Lat(KK),SS(KK+2,:))
 axis equal
       
